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Abstract 

Measurements in turbulent flows have revealed that the velocity field in nonequilib- 
rium systems exhibits ^-exponential or power law distributions in agreement with 
theoretical arguments based on nonextensive statistical mechanics. Here we consider 
Hele-Shaw flow as simulated by the Lattice Boltzmann method and find similar be- 
havior from the analysis of velocity field measurements. For the transverse velocity, 
we obtain a spatial g-Gaussian profile and a power law velocity distribution over all 
measured decades. To explain these results, we suggest theoretical arguments based 
on Darcy's law combined with the non-linear advection-diffusion equation for the 
concentration field. Power law and g-exponential distributions are the signature of 
nonequilibrium systems with long-range interactions and/or long-time correlations, 
and therefore provide insight to the mechanism of the onset of fingering processes. 
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Experiments in turbulent flows (such as Couette- Taylor flow [1]) have shown 
that data obtained from quantities extracted from velocity field measurements 
in these nonequilibrium systems exhibit power law or g-exponential velocity 
distributions which have been analyzed and interpreted with theoretical argu- 
ments based on nonextensive statistical mechanics [2]. Other physical systems 
have been shown to exhibit similar type distributions as a consequence of 
nonextensivity [3]. Here we consider fingering which occurs in the interfacial 
zone between two fluids confined between two plates with a narrow gap (Hele- 
Shaw geometry) when one of the fluids is displaced by the other fluid, the 
interfacial instability resulting from mobility differences between the fluids. 
Using a mesoscopic approach - the lattice Boltzmann method - we investi- 
gated the dynamics of spatially extended Hele-Shaw flow [4] . We are primarily 
interested in the early stage of the fingering process in order to explore the 
characteristics of the velocity field after the destabilization of the interface. 
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The Hele-Shaw geometry is intrinsically tri-dimensional, but the effective desta- 
bilization of the flow can be described in the 2-D plane, the flow being purely 
Poiseuille flow in the third dimension of the shallow layer [5]. So we can 
simulate Hele-Shaw flow by using the two-dimensional LBGK equation, the 
Bhatnagar-Gross-Krook version of the lattice Boltzmann equation where the 
full collision term is approximated by a single relaxation term [6]. One emu- 
lates 3-D flow by introducing a drag term, thereby simulating a system with 
a virtual cell gap in the third dimension: the drag enters the LBGK equation 
as a damping term [7]. 

As the simulation method has been described in detail in a previous publica- 
tion [8], we merely outline the essentials of the physics of the problem as it is 
simulated. The system consists of a 2-D rectangular box with vertical length 
L y and horizontal width L x (L x x L y = 1024 x 2048 nodes) filled with fluid 
2 which is being displaced by fluid 1 injected uniformly from the top of the 
box (the respective drag coefficients are f3\ < ($2). Fluid 1 invades the system 
through a constant drift applied along the ^-direction. The initial concentra- 
tion profile of fluid 1 is a step function which, because of mutual diffusion 
between the two miscible fluids, develops into an erfc y profile. The interfa- 
cial profile which is initially flat along the transverse x-direction, is perturbed 
with white noise to trigger the instability. Any small perturbation so induced 
is the siege of a local pressure gradient normal to the interface where from the 
instability develops as illustrated in Fig.l. 

Measurements of the time evolution of the mixing length (L mix ) of the interfa- 
cial zone (see Fig.l) show a dynamical transition from the short time diffusive 
regime to the nonlinear regime [8]. The short time behavior obeys a power 
law L m i X oc t M with /i ~ 1/2; then (as shown in Fig. 3 in [8]) a transition 
is observed from the diffusive regime with the "early time exponent" 1/2 to 
the regime with a "dynamic exponent" \i ~ 2.3, a value which was obtained 
from lattice Boltzmann simulation as well as from direct numerical simula- 
tions of Darcy's equation [8]. Here we shall restrict to the investigation of the 
transverse velocity field (v x ), that is the field transverse to the propagation 
direction of the flow, during the onset of fingering. What is meant by onset 
is the late stage of the diffusive regime with fx ~ 1/2 (see Fig. 3 in [8]) where 
the fingers have developed into a regular pattern (see Fig.l) before growing 
into the competitive stage where large fingers absorb slower smaller fingers, a 
feature characteristic of the nonlinear regime with \x ~ 2.3. 

When the planar interface destabilizes, the flow develops a wiggling concen- 
tration profile (in the x-direction) thereby producing local concentration gra- 
dients which trigger vorticity fluctuations. Pattern formation follows in the ve- 
locity field as shown in Fig. 2a, and "vortices" are being created with alternat- 
ing polarities. They are reminiscent of vortex patterns in two-dimensional tur- 
bulent flows, although here the Reynolds number has low value. However the 
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relevant control parameter, the Peclet number, is high: Pe = L x \v y \/D ~ 200 
(D is the diffusion coefficient). We performed measurements of the velocity 
field and found that the transverse component v x in each vortex shows a 
g-Gaussian profile (see Fig.2b), and that the corresponding velocity distribu- 
tion, P(v x ), obtained by computing v x over all vortices, follows a power law, 
P(v x ) (xv x - 875 (see Fig.3). 

It has been shown that processes stemming from the generalized Fokker-Planck 
equation [3,9,10] or the generalized advection-diffusion equation [11,12] can 
produce stationary solutions of the g-exponential type provided some condi- 
tions are met for the drift and dissipative functions. In the present simula- 
tions, we observe that the concentration field exhibits a g-Gaussian profile 
which suggests that the concentration C(x,t) is governed by the "porous me- 
dia equation" [13] 

d t C(x,t) + d x [K(x)C(x,t)] = d x [D(x)d x C a (x,t)], (1) 

where K(x) and D(x) are the drift and dissipation functions respectively. Note 
that the second term in (1) can be rewritten formally as 

d x [Dd x C a ] = d x [D a d x C], (2) 

where D a = aDoC a_1 . The stationary solution to Eq.(l) is a g-exponential 
of the form (see e.g. [10]) 

C(x) = C ep '*' A = C [1 - (1 - q) % \x\ x ] H , (3) 

where Co is a normalization constant, and q + a = 2. The concentration profile 
found from the simulation results has exactly this form with A = 2. 

On the other hand, the velocity field is related to the concentration through 
Darcy's law (see e.g. [14]) which can be written as 

Vx[/?(C)v] = 0, (4) 

or 

R(vx VC) + (V x v) = 0, (5) 

with R = d\n(3/dC, where (3 is the damping function which, in the fingering 
simulations, controls the drag and depends on the concentration C(x). The 
solution to Eq.(5) is a function v(C) which, through the concentration, is space 
dependent, but is analytically unknown. However, given that the concentration 
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exhibits a g-Gaussian profile (3), one can expect that the velocity field also 
assumes a g-Gaussian form. For the transverse velocity field, Fig.2b shows that 
the (/-Gaussian interpretation is in perfect agreement with the data obtained 
from the simulation results, i.e. 

v x (x) = v [1 - (1 - q)(j) q x 2 }— , (6) 



where Vq is a normalization constant. Note that the value of q in (6) is not 
the same as in (3). So the pattern of alternating vortices of the transverse 
velocity amplitude as seen in Fig. 2a can be viewed as a landscape of hills and 
wells of radial g-Gaussians distributed along the concentration profile of the 
interfacial zone. 

We next consider the velocity distribution which follows from this spatial 
profile. The g-Gaussian spatial profile (shown in Fig. 2b) is obtained experi- 
mentally by measuring the v x amplitude in a section plane orthogonal to the 
2/-axis. Now each g-Gaussian "blob" in Fig. 2a is a function of r 2 = x 2 + y 2 
and is expected to be radially symmetric. So generalizing the g-exponential in 
Eq.(6) to an arbitrary power A of the radial variable r, we have 

v(r) = v x (r)/v = e-*« rA = [1 - (1 - q)<P q r x ]^ , (7) 



We seek to establish the corresponding probability distribution function P(v) 
which, for the sake of generality, is considered in d dimensions 

oo 

P(v) = c d J r d - x \dr\ 8(v(r) - v) , (8) 
o 



where q is an integration constant. Expression (7) is inverted to give r A = 
— 4- ln q v which is used in (8) to obtain 



\r\ d x v H 



For two-dimensional systems (d = 2) with g-Gaussian spatial profile (A = 2), 
we obtain the power law distribution 

P(v) = ^v~ q . (10) 



This is exactly the power law that we find for the transverse velocity distribu- 
tion as exhibited in Fig. 3 where a fit to the simulation data gives an exponent 
q < 1 (indicating statistical importance of large values of the velocity). 
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In conclusion, we have presented an analysis and an interpretation of the 
statistical behavior of the transverse velocity field in the "onset" region of 
the fingering process just below the dynamical transition. To the best of our 
knowledge, this transition has remained unexplained analytically. The velocity 
distribution analysis developed in the present paper offers a method to inves- 
tigate both the transverse and longitudinal velocity fields, P(v x ) and P(v y ) 
respectively, below and above the dynamical transition, and should therefore 
provide insight in the nature of the transition as will be discussed elsewhere. 
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Fig. 1. Illustration of Lattice Boltzmann (BGK) simulation of viscous fingering. The 
mixing length between the two fluids is clearly identified by the shadow zone around 
the density (= 0.5) isoline (color code indicating the concentration of the invading 
fluid from red (d = 1) to blue (Ci = 0)). 
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Fig. 2. (a) Upper panel: Transverse velocity field showing landscape of g-Gaussian 
"hills and wells"; simulation data with iso- velocity contours, color code indicating 
highest positive values (red) to largest negative values (blue) of v x . (b) Lower panel: 
Transverse velocity profile (solid line) obtained from the simulation data by a section 
plane cut through the hills and wells parallel to and North- West of the main valley 
in the upper panel; analytical function Eq.(6) (dashed line): the convex and concave 
curves are connected by summation of g-Gaussians. 
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Fig. 3. Velocity distribution with power law fit : P(v) with v = (v x — (v x )) <r _1 where 
(v x ) is the mean transverse velocity and a = [(v 2 .) — (v x ) 2 \ l l 2 . Data are rescaled 
(shifted upwards in the far right branch of the upper panel) showing that all data 
follow the same power law down to values v ~ 0.05 (below which statistics become 
insufficient; here 5 points with v < 0.05 out of 1000). P(v) = A\v\~ q (Eq.(10)) with 
q = 0.8739 and A = 0.1082 for \v > 0|, and q = 0.8787 and A = 0.1054 for \v < 0|. 
Lower panel illustrates power law for 40 sets of data taken at successive times. 
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